arXiv:1503.01677vl [cond-mat.soft] 5 Mar 2015 


Insights into DNA-mediated interparticle interactions from a coarse-grained 
model 

Yajun Ding^ and Jeetain Mittal^’ 

Department of Chemical and Biomolecular Engineering, Lehigh University, Bethlehem, Pennsylvania 18015, 

USA 

(Dated: 6 March 2015) 

DNA-functionalized particles have great potential for the design of complex self-assembled materials. The 
major hurdle in realizing crystal structures from DNA-functionalized particles is expected to be kinetic bar¬ 
riers that trap the system in metastable amorphous states. Therefore, it is vital to explore the molecular 
details of particle assembly processes in order to understand the underlying mechanisms. Molecular simula¬ 
tions based on coarse-grained models can provide a convenient route to explore these details. Most of the 
currently available coarse-grained models of DNA-functionalized particles ignore key chemical and structural 
details of DNA behavior. These models therefore are limited in scope for studying experimental phenomena. 
In this paper, we present a new coarse-grained model of DNA-functionalized particles which incorporates 
some of the desired features of DNA behavior. The coarse-grained DNA model used here provides explicit 
DNA representation (at the nucleotide level) and complementary interactions between Watson-Crick base 
pairs, which lead to the formation of single-stranded hairpin and double-stranded DNA. Aggregation between 
multiple complementary strands is also prevented in our model. We study interactions between two DNA- 
functionalized particles as a function of DNA grafting density, lengths of the hybridizing and non-hybridizing 
parts of DNA, and temperature. The calculated free energies as a function of pair distance between particles 
qualitatively resemble experimental measurements of DNA-mediated pair interactions. 


I. INTRODUCTION 

It is now well recognized that self-assembly of parti¬ 
cles functionalized with biomolecules is a promising way 
to form unique nanostructured materials^’^. The dis¬ 
tinct advantage of inter-particle interactions mediated by 
DNA molecules is the specificity of Watson-Crick (WC) 
pairing^dj functionalized particles can be bridged to¬ 
gether either by direct hybridization of complementary 
single-stranded DNA (ssDNA) molecules (sticky end) or 
indirectly via “linker” DNA molecules that can bind si¬ 
multaneously to the complementary ssDNA sequences on 
two different particles^. While significant progress has 
been made in recent years, the fundamental details of 
DNA-mediated particle assembly are not very well under¬ 
stood^ Specihcally, the assembly of micron-sized par¬ 
ticles into crystalline structures is still quite challenging, 
though nanoparticles have been assembled into a wide 
variety of periodic arrangements via DNA-mediated in- 
teractions^^^. 

Computer simulations can provide a convenient 
route to explore the large parameter space of DNA- 
functionalized particles (DFPs), which can be useful for 
developing basic understanding and to test existing the¬ 
oretical design modelsUsing all-atom models to 
study the properties of DFPs is still beyond the current 
state-of-the-art computational capabilities^^^^^. Simple 
and accurate coarse-grained (CG) models are, there¬ 
fore, needed to overcome existing computational chal¬ 
lenges and are actively being developed. Starr and co- 
workers^^^^® have developed a “two-bead” DFP model 
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in which each nucleotide is represented by two force 
sites, one for the phosphate-sugar backbone and one 
for the nitrogenous base. The backbone beads are con¬ 
nected by the standard hnite extensible nonlinear elastic 
(FENE) bond potential. The beads representing the nu¬ 
cleotide bases are also connected to the backbone beads 
by the FENE potential. To model WC pairing, the 12-6 
Lennard-Jones (LJ) nonbonded potential between com¬ 
plementary bases is used. The particle itself is modeled 
as a single spherical or icosahedral core for simplicity. 
This model has extensively been applied to understand 
nanoparticle dimers, the stability of nanoparticle crys¬ 
tals, polymorphism, equilibrium clustering and dynam- 
ics^‘^“^®. Seifpour et al. also proposed a modified model 
based on this model and applied it to study the effect of 
DNA strand composition and sequence on the structure 
and thermodynamics of DFPs^®. 

Li et recently introduced a CG DFP model 

which is based on an earlier model by Travesset et . 

In this model, the DNA molecule is made up of three 
parts: single-stranded DNA, double-stranded DNA (ds- 
DNA), and a “sticky end”. Connected beads of different 
sizes are used to represent ssDNA and dsDNA. The sticky 
end is modeled by multiple beads to take into account 
the selectivity and directionality of hydrogen bonding be¬ 
tween complementary DNA bases. With this model, Li et 
al. reproduced all nine crystal structures observed exper¬ 
imentally by Macfarlane et al.^, and also proposed new 
linker sequences for future experiments. Frenkel and co- 
workers®’^®’^^’^® have developed a CG “core-blob” DFP 
model by further reducing the DNA degrees of freedom. 
In their model, the sticky end of tethered DNA is mod¬ 
eled as a blob connected to an effective spherical core. 
This model can capture several important features of 
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DNA-mediated particle assembly, from pair-particle in¬ 
teractions® to thermodynamic phase transitions^®’®®. 

When moving from the CG model with explicit DNA- 
like chains to a simplified core-blob model, the structural 
and chemical details of DNA molecules are gradually lost. 
Although it is obvious that significant coarse graining is 
needed to explore DFP system behavior (in particular 
spontaneous crystallization), certain molecular details, 
which are often ignored in these previous models, are 
important to understand DNA-mediated particle assem¬ 
bly. For example, most current models use an average 
base representation (ignoring the differences between A:T 
and G:C pair interactions), thereby neglecting the depen¬ 
dence of particle interactions on the DNA sequence. Not 
only is the difference between strengths of G:G pairing 
(stronger due to the presence of three interbase hydrogen 
bonds) compared to A:T pairing (two interbase hydro¬ 
gen bonds) important, but the actual DNA sequence can 
also affect properties of hybridized DNA molecules^®’®®. 
In addition, it has been reported that in order to capture 
surface-adsorbed or interfacial DNA structure, one needs 
to account for non-WG base pairing as well, since dsDNA 
melting and ssDNA properties will be affected by these 
additional interactions®^^®®. Base stacking interactions, 
which drive the coplanar alignment of neighboring bases, 
are also important to the overall stability of helical du¬ 
plexes and are also sequence-dependent®®. Furthermore, 
the intra-particle interactions, which are the interactions 
between DNA strands on the same particle, are often ig¬ 
nored but can be quite important for particle diffusion 
and kinetics of assembly®®. 

In this paper, we present a new GG DFP model in 
which we use a two-bead representation for each nu¬ 
cleotide of tethered DNA molecules. The DNA model, 
which is based on a previously proposed model by Dorf- 
man and co-workers®®, displays typical experimentally 
observed temperature melting behavior for dsDNA hy¬ 
bridization and ssDNA hairpin formation. The DNA 
strands are connected to spherical particle cores (made 
up of smaller point particle beads) by FENE bonds to 
simulate DEPs. We use this model to study the poten¬ 
tial of mean force (PMF) between two DEPs as a func¬ 
tion of DNA grafting density, spacer length, sticky end 
length, and temperature. The results presented in this 
paper suggest that this new DPP model can be useful 
for studying the thermodynamics and kinetics of particle 
assembly mediated by DNA hybridization. 

The paper is organized as follows: In section II, we out¬ 
line the details of the model and potential parameters. In 
section IIIA, we present data on the melting behavior of 
ssDNA hairpin and dsDNA in the absence of particles. 
Next, in section IIIB, we show potentials of mean force 
between two DFPs and discuss their connection with pre¬ 
vious experimental and theoretical data. In section IV, 
we conclude with some key observations from this work 
and discuss our future goals. 



FIG. 1. Coarse-grained DNA functionalized particle model. 
Two particles with multiple DNA molecules attached to the 
particle surface atoms are shown. The DNA sequence is com¬ 
posed of two parts, the spacer (6T: TTTTTT) and the sticky 
end as shown in the schematic. 


II. MODEL AND SIMULATION DETAILS 
A. DNA and particle models 

In order to achieve an accurate description of DNA be¬ 
havior in a GG DFP model, the DNA model itself should 
be suitable for studying the properties of ssDNA and ds¬ 
DNA. The model proposed here is expected to include 
the following key aspects of DNA behavior that may 
be important for understanding the assembly of DFPs: 
(i) attractive interactions between complementary base 
pairs A:T and G:G, (ii) attractive interactions between 
adjacent bases to capture base stacking, (iii) the ability 
to study sequence-dependent behavior, (iv) hairpin for¬ 
mation in a suitable ssDNA sequence with complemen¬ 
tary WG pairs at both ends, (v) hybridization between 
two fully or partially complementary strands, and (vi) 
no multi-strand association or aggregation in the form of 
large bundles. 

It is also desirable to include base pair dependent 
non-WG interactions®®, but is something that will re¬ 
quire extensive parameterization and testing. We plan 
to do this in future refinement of the model proposed 
here. By taking into account the requirements above, we 
have developed our GG model, which is based on an ear¬ 
lier DNA model proposed by Dorfman and co-workers®®. 
Our choice is based on a compromise between compu¬ 
tational efficiency and accurate representation of DNA 
behavior. Other DNA models are available in the lit¬ 
erature, which can be computationally more efficient or 
physically more accurate. In our model (see Figure 1), 
a two bead representation is used for a nucleotide, one 
bead for the phosphate-sugar backbone and another for 
the nitrogenous base. 

The connectivity between base and backbone beads 
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and between adjacent backbone beads is modeled by the 
finite extensible nonlinear elastic (FENE) bond poten¬ 
tial^^, which is given by 


b^bond(^ij) 
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where fcbond is the effective strength of the potential, and 
Rq is the cutoff distance at which this potential diverges. 
Following Kenward et al., we use fcbond = 30e/cr^ and 
Rq = 1.5(7 here Here, e is the characteristic energy 
scale and a is the characteristic length scale. The angle 
potential is used to provide additional bending rigidity 
to the backbone with the following functional form: 


base-base beads are modeled by the Weeks-Chandler- 
Andersen (WCA) potential as 


UwCA{rij) 
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The repulsive interactions between negatively charged 
backbone beads are modeled by the Yukawa potential 
(screened electrostatics) as^® 
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0, if Tij ^ 7 


Hangle(6') = ^^^^{cOs9 + COs9of, (2) 


where 9 the angle between three adjacent backbone 
beads, 9q = 180° is the reference bending angle, and 
fcangie = 24e is the stiffness parameter. Inter-chain and 
intra-chain stacking (st) and WC hydrogen bonding (hb) 
attractions between DNA base beads are modeled as^^ 


Ukivij) =-eukSfj ^exp (20 ^-Fs^-blj 


( 3 ) 


where k € {hb, st} and Fg sets the range of interaction; 
we use Fg = 1.5. The interactions between possible base 
pairs out of the four base alphabet (A, T, C or G) can 
be simply represented in a matrix form [6fj] as^^ 
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In Eq. (3), Uk sets the relative energy scale between 
base stacking Ust and hydrogen bonding Uhb- We use 
Ust = 2.5 and Uhb = 1-0 and keep Ugt/uhb = 2.5 as sug¬ 
gested previously by Linak et al.^^. The stacking poten¬ 
tial Ust is only applied to adjacent bases in a DNA strand 
connected to backbone beads i and f±l, which are less 
than distance 2a apart. The hydrogen bonding potential 
Ghb is applied between all possible WC inter-strand base 
pairs and intra-strand base pairs, which are less than 
distance 2a apart, excluding adjacent {i and z±l) and 
next-nearest neighbors (i and z±2). In addition, short- 
range repulsive interactions between backbone-base and 


where A is a prefactor with units of energyx distance, 
K = 2/(7 is the screening length, and rent = 3.2ct is the 
potential cutoff value. In our reduced unit model, the 
actual strength of electrostatic interactions (A = 100 ea) 
is selected so as to prevent triple- or multi-strand ag¬ 
gregation as shown in Figure 2. The undesired aggrega¬ 
tion of multiple DNA strands results from the absence of 
electrostatic repulsion between backbone beads, whereas 
well-separated dsDNA pairs form in the presence of such 
repulsion modeled as the Yukawa interaction potential 
(eq. 7). We note that aggregation of multiple DNA 
strands can also be avoided by accounting for the direc¬ 
tional nature of hydrogen bonding interactions between 
DNA base beads with a multi-body interaction poten¬ 
tial instead of the isotropic distance dependent potential 
given by eq. 3. To account for electrostatic interactions 
in a more transparent manner, we are currently devel¬ 
oping a DNA model in real units with appropriate size 
scaling for the backbone and base beads as well as dis¬ 
parate bond lengths for backbone-backbone and base- 
backbone bonds^^. No attempt is made here to convert 
the reduced simulation units to SI units as was done pre¬ 
viously by Linak et al.^^’^® for the temperature based on 
the experimental data. Without significant experimen¬ 
tal input, such a mapping will likely yield inconsistent 
results. 

The particles are modeled as rigid body hollow spheres 
(radius 5a used throughout this work), and each parti¬ 
cle is made up of 100 beads uniformly distributed on 
the surface of a sphere®®. The particle beads are un¬ 
charged and interact with DNA and other particle beads 
via the WCA interaction potential. Single strand DNA 
molecules are covalently linked to the particle surface 
beads by the FENE bond potential. We select parti¬ 
cle surface beads that are farther than a certain distance 
value, which varies with DNA grafting density, to ensure 
uniform distribution of DNA strands on the particle sur¬ 
face. Figure 1 shows a snapshot of two CG DFPs with 
several ssDNA molecules attached at random locations to 
the particle surface. In our model, DNA molecules that 
are attached to the same particle can interact via the WC 
base pair interactions, in contrast to previous CG DFP 
models that neglect intra-particle base pair interactions 
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(a) 

(b) 






FIG. 2. Simulation snapshots of four pairs of partially comple¬ 
mentary DNA sequence (see Figure 1) without (a) and with 
(b) the Yukawa interaction potential temperature below the 
DNA melting temperature. 


even between A:T and G:C base pairs. Therefore, appro¬ 
priate DNA sequences in our model can form hairpin or 
loop structures, which have been proposed as ideal can¬ 
didates for novel design schemes^®^'^^. Both inter-particle 
and intra-particle base pair interactions are naturally in¬ 
corporated in our CG DFP model. 


B. Simulation details 

To enhance equilibrium sampling, the replica exchange 
molecular dynamics (REMD) method is used for all of 
our simulations'^^. REMD simulations are carried out 
under the canonical ensemble using a Langevin thermo¬ 
stat with damping parameter r = 1 in a 

cubic simulation box with periodic boundary conditions 
in all directions. The size of the simulation box is cho¬ 
sen to be IOOct so that there are no interactions between 
molecules and their periodic images. In our REMD sim¬ 
ulations, we typically use 16 to 48 temperature replicas 
to obtain sufficient exchange probability and, therefore, 
convergence. The replica temperatures were chosen by 
trial and error such that the potential energy distribu¬ 
tions adjacent replicas overlap sufficiently to ensure ex¬ 
change between any two adjacent replicas with at least 
50% probability. The total number of timesteps that 
we performed for each simulation is of the order of 10® 
steps with a step size of At = 0.01 (e/m/cr^)”^/^. About 
10® timesteps per replica are required to obtain a rea¬ 
sonable estimate of the melting profiles of the hairpin 
and the duplex, which takes approximately 144 hours of 
processor time (for 32 replicas) on the Lehigh comput¬ 
ing server corona (corona.cc.lehigh.edu) equipped with 
AMD Opteron 8-core 6128, 2GHz processors. For DFPs, 
we need about 5 x 10® timesteps per replica, requiring ap¬ 
proximately 2304 hours of processor time (for 32 replicas) 


on the Lehigh computing server corona. 

The DNA strands (hairpin or ssDNA) or DFPs are 
placed randomly inside the simulation box. The initial 
part of the simulation is discarded as equilibration, and 
the equilibration length is decided based on observables 
such as number of hydrogen bonds and distance between 
the particles. For hairpin and dsDNA duplex simulations 
the equilibration length is less than 10^ steps, while the 
equilibration length for DFPs is about 10® steps. We use 
the block averaging scheme to estimate errors by divid¬ 
ing the production data into three sets and calculating 
the standard deviation from the mean. In case of the du¬ 
plex/hairpin, the error bars of melting curves are smaller 
than the symbol size. 

When simulating DFPs, we also tested to ensure that 
results remain unchanged with a box size of 200cr. For the 
larger box size, we find it useful to employ umbrella sam¬ 
pling simulation^® to obtain free energies as a function of 
pair distance between particles in a reasonable time. This 
also helps validate results obtained from REMD simula¬ 
tion. Eor umbrella sampling, we use harmonic umbrella 
potential with spring constant = 5.0 e/cr^ along the pair 
distance between the particles with a sufficient number of 
copies to obtain overlap between distance distributions. 

Finally, we have used two different pairs of com¬ 
plementary DNA sequences in our simulation^®’ 47; a 
12-mer, long sticky end (TTTTTTATGTATCAAGGT 
and TTTTTTAGGTTGATACAT) and a 7-mer, 
short sticky end (TTTTTTTTTTTGTGTACC and 
TTTTTTTTTTTGGTAGAG). 


III. RESULTS 

A. Melting behavior of ssDNA hairpin and dsDNA 

In order to validate the CG DNA model presented in 
section IIA, we first present the simulation data to char¬ 
acterize the melting behavior of ssDNA and dsDNA. The 
melting temperature is identified as the temperature at 
which the heat capacity (calculated from potential energy 
fluctuations) shows a maximum, the so-called calorimet¬ 
ric definition^®. We also use a structure based definition 
to define melted or unhybridized states. We define ss¬ 
DNA hairpin or dsDNA duplex configurations to be hy¬ 
bridized (see Supplemental Material (SM) Figure Sl^^) if 
at least a certain number of the possible complementary 
base pairs are bonded. Up to a certain threshold value of 
number of bonded pairs (see SM Figures S2 and S3 
the melting curves are quite similar. The melting tem¬ 
perature, identified from this structure based criteria as 
the temperature at which 50% of the states are melted, 
is similar to the calorimetric definition based on the heat 
capacity. 

Figure 3 shows the melting curves for a ssDNA hairpin: 
A 10 G 20 T 10 and a dsDNA S1S2: GGGTCATACAGTGC 
obtained from REMD simulations. The shape of these 
melting curves is similar to what is expected based on 
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FIG. 3. Thermodynamic melting behavior. The fraction of 
melted states as a function of temperature (red symbols) (and 
typical DNA configurations) and heat capacity (black sym¬ 
bols) are shown for (a) a ssDNA hairpin:AioG 2 oTio and (b) 
a pair of dsDNA S1S2:GGGTCATAGAGTGC. 


thermodynamic melting behavior of hairpins and dsDNA 
formation from experiments®*^’®^, although the widths of 
the melting transitions are much broader as found by 
Linak et al.®®’®® as well for a similar model as used here. 
Such broad transitions (unphysical) are also found in 
more detailed atomistic models of biomolecular systems 
such as proteins®^. From Figure 3 we observe that at very 
low temperatures almost all states are found to be hy¬ 
bridized and at high temperatures all states are melted. 
The transition between these two states is occurring in 
a cooperative manner, more so in the case of dsDNA. 
The temperature at which the fraction of melted states 
equals the fraction of hybridized states can be defined 
as the melting temperature. We note that the melting 
curve obtained from simulation of a single pair of du¬ 
plex in the canonical ensemble is different from the bulk 
experimental situation involving hybridization between 
many pairs of complementary strands, due to finite-size 


effects, as discussed in detail by Ouldridge et al.®®’®^ If a 
comparison between simulation data from a finite system 
and experiment is desirable, melting temperature can be 
estimated as the temperature at which 33% of states are 
melted, as opposed to 50%. We also note that in the case 
of dsDNA hybridization, the melting temperature is also 
a function of the DNA concentration or simulation box 
volume for a fixed number of DNA strands. 

Based on the results above, the CG DNA model is 
expected to capture the essential thermodynamics of ds¬ 
DNA hybridization and ssDNA properties. Specifically, 
cooperative transitions between hybridized and unhy¬ 
bridized states for dsDNA and ssDNA, as a function of 
temperature are captured. Next, we use this model to 
study interactions between a pair of DNA-functionalized 
particles. 


B. Pair interactions between DNA-functionalized particles 

Since attentive control of particle interactions is needed 
to obtain ordered structures, several previous efforts have 
been focused on understanding pair interactions between 
DNA-functionalized particles. On the experimental side, 
Crocker and co-workers measured DNA-mediated pair in¬ 
teractions between micron-sized particles.^®’®® The sepa¬ 
ration distances between particles were obtained by con¬ 
fining two DFPs in harmonic potentials in an optical 
tweezer setup. The pair interaction potential (or free 
energy as a function of distance) is then obtained by 
the Boltzmann relation, accounting for optical forces. 
These direct measurements have allowed the develop¬ 
ment of statistical physics-based theoretical models for 
understanding the underlying physical phenomena^® of 
interactions between micron-sized particles. Though ex¬ 
perimentally very challenging, it is desirable to obtain 
such information for interactions between nanoparticles 
as well. 

Estimates of pairwise interactions between nanopar¬ 
ticles have been obtained by simulation with the help 
of CG DFP models (see section 1)56.57 xj^ese previous 
studies consider pair interactions between DFPs for very 
low DNA grafting densities of 4 to 6 DNA molecules per 
particle. Frenkel and co-workers®®’®® have developed nu¬ 
merical schemes to extract pair interactions using their 
core-blob DFP model with a combination of statistical 
mechanics theory and Monte Carlo simulations. Our new 
DFP model presented here is meant to act as a bridge be¬ 
tween the core-blob type simplified description (which is 
computationally tractable even for phase behavior cal¬ 
culations) and experiment by including a more detailed 
description of DNA degrees of freedom (see section IIA). 
Ultimately, it may turn out that some of these details 
are not necessary to capture essential details of DFP as¬ 
sembly thermodynamics, including phase diagrams. But 
we anticipate that some of these additional details in our 
model may be important for understanding issues related 
to particle dynamics and the kinetics of assembly. 
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We simulate two DFPs grafted with either the 12-mer 
long sticky end or the 7-mer short sticky end as described 
previously. The potential of mean force (PMF) or free 
energy as a function of pair distance between two DFPs, 
with respect to unhybridized particles, is calculated as 

AF/(r) = -kBT\D.[g{r)] + c, (8) 

where g{r) is the pair correlation function (normalized 
with respect to an ideal gas), fcs is the Boltzmann con¬ 
stant, c is an additive constant, and I is used to identify 
particle types in a pair, i.e., I = {A-A, B-B, A-B}. In ad¬ 
dition to binning data collected at a single temperature 
in an REMD simulation, we also use the weighted his¬ 
togram analysis method (WHAM)®®’®° to combine data 
from all of the temperatures to obtain PMFs. 




20 30 40 50 


Pair distance [a] 

FIG. 4. Potential of mean force (PMF) calculated for a graft¬ 
ing density of 16 DNA molecules per particle with long sticky 
end for pairs A-B (a), A-A (b, top panel), and B-B (b, bottom 
panel). The data from simple binning at a single temperature 
shown by empty symbols, and solid lines are estimates from 
the weighted histogram analysis method. Symbol and line 
colors are used to distinguish between various temperatures 
as labeled in panel a. 


1. Effect of temperature 

Figure 4 shows a representative set of PMFs for dif¬ 
ferent temperatures (see legend) between three possible 
pairs of particles (A-A, B-B, and A-B). The symbols rep¬ 
resent data from the REMD simulations sampled at a 
single temperature, and the lines are the WHAM re¬ 
construction based on the data from all temperatures. 
The agreement between the two estimates clearly indi¬ 
cates that the simulation times are long enough to ob¬ 
tain PMFs. The interactions between like pairs (A-A 
and B-B) are purely repulsive at all of the temperatures, 
as shown in Fig. 4b, as these ssDNA sequences were de¬ 
signed to have minimal self-interactions. Moreover, the 
repulsive interactions between like pairs are independent 
of the temperature, except at very low pair distances. 

The PMFs between unlike particle pairs (A-B), which 
are coated with complementary ssDNA with a long sticky 
end, are shown in Fig. 4a. The free energy as a func¬ 
tion of pair distance shows a distinct minimum at low 
temperatures, which corresponds to the bound configu¬ 
ration between particles A and B (stabilized by DNA 
hybridization). Specifically, the general shape of the pair 
potential and its dependence on temperature (compare 
Fig. 4a with Fig. IE from Ref. 46) resemble the exper¬ 
imentally measured effective pair potentials^®. For two 
particle simulations in a large enough box (approaching 
zero density limit), the minimum in free energy is an im¬ 
portant parameter to characterize interactions between 
particles as a function of system parameters such as tem¬ 
perature, DNA grafting density, and length of sticky or 
spacer end®^. 

Figure 5 shows the minimum interaction free energy 
as a function of temperature, which we find to vary in 
a highly non-linear manner. Previously, Dreyfus et al.^^ 
also predicted a similar non-linear dependence of the min¬ 
imum interaction free energy on temperature with the 
help of a theoretical model. In their model, the mini¬ 
mum interaction free energy depends on the number of 
bonds formed between the two particles as well as the 
entropy associated with the selection of various combi¬ 
nations of strands or bonds on the two particles; the 
sharp decrease in the minimum interaction free energy 
as a function of decreasing temperature is most likely 
due to enhanced hydrogen bonding (due to the increase 
in the number of hydrogen bonds as the temperature de¬ 
creases and the increased stability of a hydrogen bond 
with decreasing temperature), which is further facilitated 
by shorter distances at which the minimum occurs. How¬ 
ever, in their model DNA is simplified as a rigid rod with 
a sticky point at the end (to represent complementary 
ssDNA), and the two particles are approximated as flat 
plates. In our model, ssDNA strands are relatively flexi¬ 
ble, and the length of the sticky end is comparable to the 
length of the spacer. Given these differences between our 
simulation model and the theoretical model of Dreyfus 
et al., the non-linear dependence of the minimum inter¬ 
action free energy on temperature is expected to be a 
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FIG. 5. The minimum interaction free energy is shown as a 
function of temperature corresponding to system parameters 
used in Fig. 4a. The inset shows a zoomed-in view of the main 
plot. 


general feature of DNA-mediated particle interactions. 
At low temperatures, the minimum interaction free en¬ 
ergy can be quite significant compared to the thermal 
energy, thereby leading to irreversible particle binding. 
Appearance of such irreversible particle binding at tem¬ 
peratures not too far from the desired temperature range 
(weak binding regime) due to non-linearity is likely to be 
an inherent hindrance to particle rearrangement, and, 
therefore, crystallization in DNA-mediated particle as¬ 
sembly. 


2. Effect of DNA grafting density 

In experiments, DNA grafting density on the particle 
surface can be controlled by synthesis methods®^’®^ as 
well as by adding non-hybridizing ssDNA in the solution 
buffer®"^. By controlling the fraction of sticky ends that 
can hybridize with complementary ssDNA on other parti¬ 
cles, it was shown that the particle dissociation behavior 
(melting temperature and sharpness) can be controlled 
quantitatively®"*^’®^. It is still less well understood how 
this particular parameter can affect the pair interaction 
potential and the associated assembly mechanism. Here, 
we study the first question by calculating PMFs as a 
function of DNA grafting density ranging from 1 DNA 
molecule per particle to 32 DNA molecules per particle. 

Figure 6 shows the minimum interaction free energy 
as a function of temperature for varying DNA grafting 
densities for the case of a short sticky end. For a given 
temperature, the minimum interaction free energy be¬ 
comes lower with increasing DNA grafting densities as 
a greater number of complementary DNA molecules can 
hybridize between a pair of particles. According to Drey¬ 
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FIG. 6. The minimum interaction free energy is shown as 
a function of temperature for different grafting densities (see 
legend) for the short sticky end. 

fus et al.®*^, the non-linear dependence of the minimum 
interaction free energy can be categorized into three dif¬ 
ferent regimes. The weak-binding regime, in which the 
minimum interaction free energy is comparable to the 
thermal energy, is most relevant for successful DNA- 
mediated particle assembly in laboratory experiments. 
In this regime, the minimum interaction free energy is 
expected to be proportional to the average number of 
hydrogen bonds formed ®*^. In Fig. 6, we limit the range 
of the minimum interaction free energy from —IQksT 
to —ksT, which approximately corresponds to the so- 
called weak-binding regime. We find that even in this 
regime, the variation in minimum interaction free energy 
as a function of temperature is quite non-linear. Similar 
behavior in the weak binding regime was also observed 
previously by Mirjam et al.^. 

To identify changes in pair PMF as a function of DNA 
grafting density, we use two special cases - constant tem¬ 
perature and constant minimum interaction free energy. 
Figure 7a shows PMFs between unlike particles (A-B) for 
different grafting densities, but the minimum interaction 
free energy is kept constant at -Ak-^T by changing the 
temperature (see legend). We find that the change in 
temperature is related logarithmically to the DNA graft¬ 
ing density in order to obtain a constant minimum inter¬ 
action free energy (data not shown). Though not entirely 
unexpected, the DNA grafting density significantly alters 
the shape of the pair interaction potential or PMF. The 
minimum in free energy is quite broad for lower graft¬ 
ing densities, but the attractive well becomes narrower 
with increasing DNA grafting density. This narrowing of 
the attractive well is most likely due to stronger repul¬ 
sion caused by overlap between grafted DNA molecules 
at small pair distances. Although some previous work 
has been done relating shape of the pair potential with 
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particle dynamics®®, there is no general framework to de¬ 
termine which of these potential shapes will be most de¬ 
sirable for ordered particle assembly. Our hypothesis is 
that intermediate DNA grafting densities will yield an 
optimum in the width of the attractive well so that pair 
distance can change as needed for structural rearrange¬ 
ments, while the resulting ordered structure is still stable 
without significant fluctuations in lattice parameters. 

In Figure 7b, we compare PMFs at the same temper¬ 
ature for varying DNA grafting densities. The minimum 
interaction free energy decreases with increasing DNA 
grafting density as expected due to the higher number 
of complementary strands that can hybridize. However, 
the maximum number of hydrogen bonds, defined as the 
number of hydrogen bonds formed between the two par¬ 
ticles at the lowest temperature of our simulation, does 
not linearly increase with the number of grafted DNA 
chains (see SM Figure 84"*®). This is most likely due to 
the excluded volume (repulsive) interactions between the 
DNA molecules on the same particle, thereby limiting 
the number of strands that can hybridize with increas¬ 
ing DNA grafting density. Similar to Figure 7a, we also 
observe narrowing of the attractive well in addition to 
deepening at higher DNA grafting densities, again due 
to enhanced repulsion between DNA molecules at small 
pair distances. 

Finally, we note that the current model, due to its ap¬ 
proximate nature, can only provide qualitative changes 
in the potential function, and therefore the observed 
changes in the potential width should be interpreted ac¬ 
cordingly. 


3. Effect of spacer length 


Experimentally, Jin et showed that an increase 
in the spacer length (the non-hybridizing part of the ss- 
DNA) increases the melting temperature of DNA grafted 
on nanoparticles. Later, Sun et o7®® also observed a simi¬ 
lar increase in DNA melting temperature with increasing 
spacer length. From molecular simulations of the CG 
DFP model, Fernando et al}^ found that longer spacer 
length can actually destabilize a preformed crystal lat¬ 
tice of DFPs at lower temperatures. Although there are 
several studies on the relationship between spacer length 
and DNA melting temperature, relatively little informa¬ 
tion is available on the dependence of the pair interaction 
potential between DFPs on spacer length. 

Figure 8 shows PMFs between a pair of unlike DFPs 
(A-B) when each particle is grafted with 16 complemen¬ 
tary ssDNA molecules with the long sticky end at the 
same temperature. The minimum interaction free en¬ 
ergy decreases slightly (by about 1 k^T) with increasing 
spacer lengths from 0 to 8 bases. This suggests that in¬ 
creasing the spacer length may help in the formation of 
more hydrogen bonds between DFPs, but other param¬ 
eters such as DNA grafting density and particle curva¬ 
ture may be more important determinants of enhanced 




Pair distance [a] 

FIG. 7. (a) Potential of mean force (PMF) for different DNA 
grafting densities for the short sticky end with the same at¬ 
traction free energy of —JfcsT. (b) PMF for different grafting 
densities for the short sticky end at the same temperature of 
0.23. 


hydrogen bond formation. More importantly, the pair 
distance at which the minimum in the PMF occurs shifts 
to a higher value with increasing spacer length. These 
results are consistent with experimental observations, as 
the decrease in the minimum interaction free energy will 
lead to a higher melting temperature and the increase 
in preferred pair distance will lead to a longer transla¬ 
tional lattice parameter^. As the width of the PMF also 
increases with increasing spacer length, one may expect 
destabilization of certain crystal lattices as observed by 
Fernando et in their simulations. 


4. Effect of sticky end length 

In addition to varying the DNA grafting density, one 
can also alter the minimum interaction free energy be¬ 
tween DFPs by changing the number of complementary 
base pairs (or type of base pairs formed, A:T versus 
G:C^®) of the sticky end. Crocker et al.^^ found that re¬ 
ducing the sticky end length leads to a decrease in melt- 
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FIG. 8. Effect of spacer length on the potential of mean force 
for long sticky end with 16 DNA molecules per particle at 
temperature 0.224. 



Pair distance [a] 


FIG. 9. Gomparison of potential of mean force between short 
and long sticky ends at three different temperatures, ss: short 
sticky end; Is: long sticky end. 


ing temperature. Similar conclusions were also drawn 
by others using separate linker-mediated interactions be¬ 
tween DFPs^°. Figure 9 compares PMFs obtained for a 
pair of DFPs with short and long sticky ends for three 
temperatures. Both of the DNA sequences (long and 
short sticky ends) have a total of 18 bases, out of which 
7 and 12 are part of the sticky end for short and long 
sticky ends, respectively. Intuitively, varying the length 
of the sticky end (while keeping the total DNA length 
constant) can have two effects on the shape of the PMF 
- change in the depth and change in the width of the at¬ 
tractive potential well. As shown in Fig. 9, the minimum 
interaction free energy for the short sticky end is higher 


than for the long sticky end at the same temperature, 
as expected. The minimum in free energy for the short 
sticky end is also shifted to longer pair distances due to 
shorter overlap between shorter sticky ends as compared 
to long sticky ends. Interestingly, the PMF well in the 
case of the short sticky end is wider as compared to that 
of the long sticky end. Current CG DFP models, which 
assume a point-like sticky end, cannot capture the effect 
of sticky end length on the shape of the PMF. This is 
not expected to be limiting for small variations in sticky 
end length but may be important for large changes in 
the sticky end length to model associated changes in the 
particle assembly behavior. 


IV. CONCLUSIONS 

In this paper, we present a new coarse-grained model 
for studying the behavior of DNA-functionalized parti¬ 
cles. The coarse-grained DNA model used here provides 
explicit DNA representation (at the nucleotide level) and 
complementary interactions between Watson-Crick base 
pairs, which lead to the formation of ssDNA hairpins 
and dsDNA. Aggregation between multiple complemen¬ 
tary strands is prevented. We use this model to calcu¬ 
late the dependence of the free energy on the distance 
between a pair of DNA-coated particles as a function of 
temperature, DNA grafting density, and lengths of sticky 
end and spacer DNA. The change in the minimum inter¬ 
action free energy as a function of system temperature is 
found to be non-linear even in the weak-binding regime. 
Our results of particle pair potentials as a function of 
system parameters can also be useful for future design of 
crystal lattices based on DNA-mediated particle assem¬ 
bly. For example, Macfarlane et al7, suggests that the 
lattice spacing can be adjusted by varying the hydrody¬ 
namic size ratio between two types of particles grafted 
with complementary DNA molecules. Based on our ef¬ 
fective pair potentials, one should be able to modify the 
lattice spacing with the same hydrodynamic particle size 
ratio by varying the sticky end length while keeping the 
total DNA length the same. Moreover, we find that sev¬ 
eral system parameters (such as DNA grafting density 
and sticky end length) can be used to change the width 
of the attractive well. We hypothesize that shallower and 
broader attractive potentials may be preferable for parti¬ 
cle assembly in ordered structures due to dynamic DNA 
hybridization^^, but may also lead to lower lattice sta¬ 
bility. In the future, we plan to address these questions 
with multiparticle assembly simulations using the new 
coarse-grained model. 
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